Loading libraries

library(dada2); packageVersion("dada2")
## [1] '1.18.0'
library(ggplot2); packageVersion("ggplot2")
## [1] '3.3.3'
library(phyloseq); packageVersion("phyloseq")
## [1] '1.34.0'
library(data.table); packageVersion("data.table")
## [1] '1.13.6'
library(dplyr)
library(plyr)
library(magrittr)
library(RColorBrewer)
library(picante)
library(igraph)
library(cowplot)
library(plotly)
library(lemon)
library(GGally)
library(randomcoloR)
library(knitr)
library(kableExtra)
library(ggpmisc)
library(ggpubr)
library(decontam)
library(stringr)

Importing data

SVt = readRDS("../seqtab_final.rds")
# SVt_cnn = collapseNoMismatch(SVt, minOverlap = 20,
#                               identicalOnly = FALSE, vec = TRUE, verbose = TRUE)
# saveRDS(SVt_cnn, "../seqtab_final_cnn.rds")
# SVt = readRDS("../seqtab_final_cnn.rds")
SVt = SVt[!rownames(SVt) %in% "Undetermined", ]

Tax = readRDS("../tax_final.rds")

metadata<-read.csv("../../220319_for_graph.csv", header=TRUE)
# metadata[, 10:47]  <- as.numeric(metadata[, 10:47])
metadata$Snowpack_Layer = gsub("Glacial Ice", "GL ICE", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Sup. Ice", "SUP ICE", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Top", "TOP", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Mid", "MID", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Basal", "BASAL", metadata$Snowpack_Layer)

order<-rownames(SVt)
metadata$Sample_ID=as.character(metadata$Sample_ID)
metadata_ord<-metadata[match(order, metadata$Sample_ID),]
rownames(metadata_ord) <- metadata_ord$Sample_ID
metadata_ord$Sample_Name <- factor(metadata$Sample_Name, levels=unique(metadata$Sample_Name))
ps <- phyloseq(otu_table(SVt, taxa_are_rows=FALSE),
               sample_data(metadata_ord),
               tax_table(Tax))
ps.b <- prune_taxa(taxa_sums(ps) > 0, ps)
ps.b = subset_taxa(ps.b, Kingdom %in% c("Archaea", "Bacteria"))
ps.b = subset_samples(ps.b, Fox_Aspect != "N_NE_D")

Testing batch effects

#remove samples with 0 reads
ps.batch = prune_samples(sample_sums(ps.b)>=1, ps.b)

readsumsdf_samples<-data.frame(nreads = sample_sums(ps.batch),
                               samples = rownames(as.data.frame(sample_sums(ps.batch))),
                               batch = sample_data(ps.batch)$Batch.no.)
readsumsdf_samples_f = readsumsdf_samples[- grep("Control", readsumsdf_samples$batch),]

my_comparisons=list(c("1","2"),c("1","3"),c("1","4"),c("1","5"),
                    c("2","3"), c("2","4"), c("2","5"),
                    c("3","4"), c("3","5"),
                    c("4","5"))

ggplot(readsumsdf_samples_f, aes(x=batch, y=nreads,
                                 color=batch, fill=batch)) +
  geom_violin(alpha=0.3) +
  geom_jitter(shape=16, position = "jitter") +
  stat_compare_means(comparisons = my_comparisons,
                     method = "t.test",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Batch Number") +
  guides(fill=FALSE) +
  xlab("Extraction batch") +
  ylab("Number of reads obtained") +
  theme(axis.text.x = element_text(vjust = 1,
                                   hjust = 1,
                                   size = 12),
        strip.text = element_text(size = 12, margin = margin()))

1 - Library sizes

sort(sample_sums(ps.b))
##     C05-cDNA-29     D12-cDNA-MM     C10-cDNA-34 A02-cDNA-MM-H20     C08-cDNA-32 
##               0               0               5               7              25 
##     C07-cDNA-31     B07-cDNA-17 E01-cDNA-MM-H20     B11-cDNA-22     B09-cDNA-19 
##              27              43              55             123             187 
##     D01-cDNA-37     A01-cDNA-MM      A11-cDNA-9     D10-cDNA-46     B10-cDNA-21 
##             262             290            2320            3844            9170 
##      A10-cDNA-8     B05-cDNA-15       B08-DNA-4       B11-DNA-7     B08-cDNA-18 
##           15595           15781           16781           29350           29808 
##     C03-cDNA-27      A08-cDNA-6      A06-cDNA-4      A09-cDNA-7      A07-cDNA-5 
##           31806           33041           37485           40604           44392 
##       A11-DNA-9      A05-cDNA-3       A08-DNA-6     B02-cDNA-12       D08-DNA-8 
##           50662           52562           53872           57398           57493 
##      E02-DNA-MM     C02-cDNA-26     D07-cDNA-43       D06-DNA-6       D02-DNA-2 
##           58537           59570           59805           60653           61230 
##      A01-DNA-MM       B12-DNA-8     B12-cDNA-23      A12-DNA-10       C07-DNA-5 
##           61807           62146           62235           63891           64148 
##     B04-cDNA-14       A07-DNA-5  A02-DNA-MM-H20      A04-cDNA-2     B03-cDNA-13 
##           68029           69142           69499           70613           76895 
##       D09-DNA-9      B02-DNA-12       A05-DNA-3     C04-cDNA-28       C05-DNA-3 
##           80028           83134           84732           86983           87568 
##       C04-DNA-2     D06-cDNA-42       A09-DNA-7     B01-cDNA-11      B04-DNA-14 
##           96798          101335          108915          110464          113222 
##    C03-DNA-B2-1      B03-DNA-13       C08-DNA-6     A12-cDNA-10     B06-cDNA-16 
##          114048          118548          122463          123046          123687 
##     C12-cDNA-36     D02-cDNA-38       D03-DNA-3       C06-DNA-4  E03-DNA-MM-H20 
##          123719          137347          138813          139808          140502 
##       D07-DNA-7      C12-DNA-10       D04-DNA-4       B09-DNA-5     D08-cDNA-44 
##          140668          148596          151009          155211          158902 
##       A06-DNA-4       C01-DNA-9       B10-DNA-6       A10-DNA-8       C09-DNA-7 
##          166098          169126          170463          170498          173577 
##      B01-DNA-11       E01-DNA-4     C06-cDNA-30     C09-cDNA-33     D03-cDNA-39 
##          174897          175874          188016          190708          192008 
##      A03-cDNA-1       A03-DNA-1       B06-DNA-2     C11-cDNA-35       A04-DNA-2 
##          193167          196436          198649          198996          200229 
##       B07-DNA-3     D04-cDNA-40       C10-DNA-8     D11-cDNA-47    D01-DNA-B4-1 
##          204577          204668          212618          220186          238199 
##     D09-cDNA-45       D11-DNA-2       D05-DNA-5    D10-DNA-B5-1     C01-cDNA-25 
##          240324          252560          267031          280585          293592 
##     D05-cDNA-41       D12-DNA-3 
##          302761          312011
df <- as.data.frame(sample_data(ps.b)) # Put sample_data into a ggplot-friendly data.frame
df$LibrarySize <- sample_sums(ps.b)
df <- df[order(df$LibrarySize),]
df$Index <- seq(nrow(df))
ggplot(data=df, aes(x=Index, y=LibrarySize, color=Sample_type)) + geom_point()

2 - Evaluating dataset contamination with decontam

#remove samples with 0 reads
ps.b = prune_samples(sample_sums(ps.b)>=1, ps.b)

sample_data(ps.b)$is.neg <- sample_data(ps.b)$Sample_type == "CONTROL"
contamdf0.5.prev <- isContaminant(ps.b, method="prevalence", threshold = 0.1,
                               neg="is.neg")

ps.final.noncontam <- prune_taxa(!contamdf0.5.prev$contaminant, ps.b)
ps.final.noncontam.noconts = subset_samples(ps.final.noncontam, Sample_type != "CONTROL")

ps.final.noncontam.noconts.nodoubs = subset_samples(ps.final.noncontam.noconts,
    Sample_ID != "D02-DNA-2" & Sample_ID != "B11-DNA-7" & Sample_ID != "C10-cDNA-34" & Sample_ID != "A08-cDNA-6")

ps.final.noncontam.noconts.nodoubs = prune_samples(sample_sums(ps.final.noncontam.noconts.nodoubs)>=1000,
                                                   ps.final.noncontam.noconts.nodoubs)

ps.final <- prune_taxa(taxa_sums(ps.final.noncontam.noconts.nodoubs) > 0, ps.final.noncontam.noconts.nodoubs)
ps.final = prune_samples(sample_sums(ps.final)>=1, ps.final)

ps.final.r <-  transform_sample_counts(ps.final, function(x) x / sum(x))
ps.final.r.f = filter_taxa(ps.final.r, function(x) sum(x) > .001, TRUE)
ps.final.r
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 35534 taxa and 65 samples ]
## sample_data() Sample Data:       [ 65 samples by 47 sample variables ]
## tax_table()   Taxonomy Table:    [ 35534 taxa by 6 taxonomic ranks ]
ps.final.r.f
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2886 taxa and 65 samples ]
## sample_data() Sample Data:       [ 65 samples by 47 sample variables ]
## tax_table()   Taxonomy Table:    [ 2886 taxa by 6 taxonomic ranks ]
ps.final
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 35534 taxa and 65 samples ]
## sample_data() Sample Data:       [ 65 samples by 47 sample variables ]
## tax_table()   Taxonomy Table:    [ 35534 taxa by 6 taxonomic ranks ]
sum(rowSums(otu_table(ps.final)))
## [1] 6750655
gl_ice_spls = as.character(unique(sample_data(ps.final)[str_detect(sample_data(ps.final)$Sample_Name, "Gl ice"),]$Sample_Name))

Setting order for all plots

sample_data(ps.final.r)$Month = factor(sample_data(ps.final.r)$Month, 
                                 levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r)$Fox_Aspect = factor(sample_data(ps.final.r)$Fox_Aspect, 
                                 levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r)$Snowpack_Layer = factor(sample_data(ps.final.r)$Snowpack_Layer, 
                                 levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))

sample_data(ps.final)$Month = factor(sample_data(ps.final)$Month, 
                                 levels = c("April","June","Early July", "Late July"))
sample_data(ps.final)$Fox_Aspect = factor(sample_data(ps.final)$Fox_Aspect, 
                                 levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final)$Snowpack_Layer = factor(sample_data(ps.final)$Snowpack_Layer, 
                                 levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))

sample_data(ps.final.r.f)$Month = factor(sample_data(ps.final.r.f)$Month, 
                                 levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r.f)$Fox_Aspect = factor(sample_data(ps.final.r.f)$Fox_Aspect, 
                                 levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r.f)$Snowpack_Layer = factor(sample_data(ps.final.r.f)$Snowpack_Layer, 
                                 levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))

3a - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Fox_Aspect, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed

ps.final.r.glom <- tax_glom(subset_samples(ps.final.r, 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)), taxrank = 'Phylum', NArm = FALSE)
ps.final.r.glom.psdf <- data.table(psmelt(ps.final.r.glom))
ps.final.r.glom.psdf$Phylum <- as.character(ps.final.r.glom.psdf$Phylum)

ps.final.r.glom.psdf[, mean := mean(Abundance, na.rm = TRUE), 
                 by = "Phylum"]
ps.final.r.glom.psdf[(mean <= 0.01), Phylum := "Other"]
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Bacteria"] <- "Unc. Bacteria"
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Archaea"] <- "Unc. Archaea"
#Removing Proteobacteria from previous table
ps.final.r.glom.psdf.noProteo<-ps.final.r.glom.psdf[!grepl("Proteobacteria", ps.final.r.glom.psdf$Phylum),]

#New table, with Proteobacterial classes only
ps.final.r.ProtCl = subset_taxa(ps.final.r, Phylum == "Proteobacteria")
ps.final.r.ProtClglom <- tax_glom(ps.final.r.ProtCl, taxrank = 'Class', NArm = FALSE)
ps.final.r.ProtClglom.psdf <- data.table(psmelt(ps.final.r.ProtClglom))
ps.final.r.ProtClglom.psdf$Class <- as.character(ps.final.r.ProtClglom.psdf$Class)

ps.final.r.ProtClglom.psdf[, mean := mean(Abundance, na.rm = TRUE), 
                       by = "Class"]
ps.final.r.ProtClglom.psdf[(mean <= 0.001), Class := "Other Proteobacteria"]
ps.final.r.ProtClglom.psdf.noP <- subset(ps.final.r.ProtClglom.psdf, select = -Phylum)

names(ps.final.r.glom.psdf.noProteo)[names(ps.final.r.glom.psdf.noProteo) == 'Phylum'] <- 'Taxa'
names(ps.final.r.ProtClglom.psdf.noP)[names(ps.final.r.ProtClglom.psdf.noP) == 'Class'] <- 'Taxa'
# ps.final.r.glom.psdf.noProteo = subset(ps.final.r.glom.psdf.noProteo, select=-c(Class))
ps.final.r.ProteoClAllPhy <- rbind(ps.final.r.glom.psdf.noProteo, ps.final.r.ProtClglom.psdf.noP)

tol18rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#771122", "#AA4455")

ps.final.r.ProteoClAllPhy$Taxa = factor(ps.final.r.ProteoClAllPhy$Taxa, 
              levels = c("Alphaproteobacteria", "Deltaproteobacteria",
                         "Gammaproteobacteria",
                         "Other Proteobacteria",
                         "Actinobacteria",
                         "Acidobacteria","Bacteroidetes",
                         "Chloroflexi","Cyanobacteria",
                         "Firmicutes","Gemmatimonadetes","Other"))
ggplot(ps.final.r.ProteoClAllPhy, 
       aes(x = Month, y = Abundance, fill = Taxa)) + 
  facet_grid(DNA_or_cDNA~Fox_Aspect,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill") +
  scale_fill_manual(values = tol18rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 12),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 12)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
  ylab("Relative Abundance\n")

3b - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Snowpack_Layer

ggplot(ps.final.r.ProteoClAllPhy, 
                                   aes(x = Month, y = Abundance, fill = Taxa)) + 
  facet_grid(DNA_or_cDNA~Snowpack_Layer,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill") +
  scale_fill_manual(values = tol18rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 11),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 12)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
  ylab("Relative Abundance\n")

3c - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Month

tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#AA4455")

ggplot(ps.final.r.ProteoClAllPhy,
                                   aes(x = Sample_Name, y = Abundance, fill = Taxa)) +
  facet_grid(DNA_or_cDNA~Month,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) +
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

4 - Top SVs and their sequences

Top20SVs_names = names(sort(taxa_sums(ps.final.r), TRUE)[1:20])
Top20SVs_names_df=as.data.frame(as(tax_table(prune_taxa(Top20SVs_names, ps.final.r)), "matrix")[,2:6])
rownames(Top20SVs_names_df) = c()
Top20SVs_names_df
##            Phylum               Class                 Order            Family
## 1  Proteobacteria Alphaproteobacteria           Rhizobiales  Beijerinckiaceae
## 2  Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 3  Proteobacteria Alphaproteobacteria           Rhizobiales         Labraceae
## 4  Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 5  Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 6  Proteobacteria Alphaproteobacteria           Rhizobiales Xanthobacteraceae
## 7  Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 8  Proteobacteria Alphaproteobacteria           Rhizobiales  Beijerinckiaceae
## 9  Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 10 Proteobacteria Alphaproteobacteria       Acetobacterales  Acetobacteraceae
## 11 Actinobacteria      Actinobacteria         Micrococcales Microbacteriaceae
## 12 Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 13 Proteobacteria Alphaproteobacteria      Sphingomonadales Sphingomonadaceae
## 14 Proteobacteria Gammaproteobacteria Betaproteobacteriales  Burkholderiaceae
## 15 Actinobacteria      Actinobacteria         Micrococcales    Micrococcaceae
## 16 Actinobacteria      Actinobacteria     Corynebacteriales      Nocardiaceae
## 17 Actinobacteria      Actinobacteria            Frankiales              <NA>
## 18 Actinobacteria      Actinobacteria         Micrococcales Cellulomonadaceae
## 19 Actinobacteria      Actinobacteria            Frankiales              <NA>
## 20 Actinobacteria      Actinobacteria            Frankiales              <NA>
##                                         Genus
## 1                                       Bosea
## 2                                 Cupriavidus
## 3                                      Labrys
## 4                                Sphingomonas
## 5  Burkholderia-Caballeronia-Paraburkholderia
## 6                              Bradyrhizobium
## 7                                Sphingomonas
## 8                            Methylobacterium
## 9                                Sphingomonas
## 10                                       <NA>
## 11                            Salinibacterium
## 12                                   Massilia
## 13                               Sphingomonas
## 14                                    Delftia
## 15                               Arthrobacter
## 16                                       <NA>
## 17                                       <NA>
## 18                                       <NA>
## 19                                       <NA>
## 20                                       <NA>
rownames(as.data.frame(as(tax_table(prune_taxa(Top20SVs_names, ps.final.r)), "matrix")[,2:6]))
##  [1] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTGTCCGGGAAGATAATGACTGTACCGGAAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGCGCGTAGGCGGACTCTTAAGTCGGGGGTGAAAGCCCAGGGCTCAACCCTGGAATTGCCTTCGATACTGGGAGTCTTGAGTTCGGAAGAGGTTGGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCGGTGGCGAAGGCGGCCAACTGGTCCGAAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"                         
##  [2] "TGGGGAATTTTGGACAATGGGGGCAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAATCGCGCTGGTTAATACCTGGCGTGGATGACGGTACCGGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTTGTAAGACAGGCGTGAAATCCCCGGGCTTAACCTGGGAATTGCGCTTGTGACTGCAAGGCTAGAGTGCGTCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCCCCTGGGACGTGACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
##  [3] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGACGGCCTTAGGGTTGTAAAGCTCTTTTAACAGGGACGATAATGACGGTACCTGTAGAATAAGCCCCGGCAAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGATTGTTAAGTCGGGGGTGAAATCCTGAGGCTCAACCTCAGAACTGCCTTCGATACTGGCGATCTTGAGTTCGGAAGAGGTTGGTGGAACAGCTAGTGTAGAGGTGAAATTCGTAGATATTAGCTAGAACACCAGTGGCGAAGGCGGCCAACTGGTCCGATACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
##  [4] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCTGGTGCTCAACACCAGAACTGCCTTTTAGACTGCATCGCTTGAATCCAGGAGAGGTGAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
##  [5] "TGGGGAATTTTGGACAATGGGGGAAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAAACCTTTGCACTAATACTGTGAGGGGATGACGGTACCGGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTCGTTAAGACAGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTTGTGACTGGCGAGCTAGAGTATGGCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCCCCTGGGCCAATACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
##  [6] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTTGTGCGGGAAGATAATGACGGTACCGCAAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGTGCGTAGGCGGGTCTTTAAGTCAGGGGTGAAATCCTGGAGCTCAACTCCAGAACTGCCTTTGATACTGAAGATCTTGAGTTCGGGAGAGGTGAGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCAGTGGCGAAGGCGGCTCACTGGCCCGATACTGACGCTGAGGCACGAAAGCGTGGGGAGCAAACA"                         
##  [7] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCCAGAGCTCAACTCTGGAATTGCCTTTTAGACTGCATCGCTTGAATCATGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACATGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
##  [8] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTATCCGGGACGATAATGACGGTACCGGAGGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGCGCGTAGGCGGCGTTTTAAGTCGGGGGTGAAAGCCTGTGGCTCAACCACAGAATGGCCTTCGATACTGGGACGCTTGAGTATGGTAGAGGTTGGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCGGTGGCGAAGGCGGCCAACTGGACCATCACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"                         
##  [9] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTTACCCGGGAAGATAATGACTGTACCGGGAGAATAAGCCCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTCAGAGGTGAAAGCCTGGAGCTCAACTCCAGAACTGCCTTTGAGACTGCATCGCTTGAATCCAGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
## [10] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGTCTTCGGATTGTAAAGCACTTTCGGCGGGGACGATGATGACGGTACCCGCAGAAGAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATGACTGGGCGTAAAGGGCGCGTAGGCGGTTGCATAAGTTAGATGTGAAATTCCCGGGCTTAACCTGGGGACTGCATTTGATACTATGTGGCTTGAGTATGGAAGAGGGTCGTGGAATTCCCAGTGTAGAGGTGAAATTCGTAGATATTGGGAAGAACACCGGTGGCGAAGGCGGCGACCTGGTCCATAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"                         
## [11] "TGGGGAATATTGCACAATGGGCGCAAGCCTGATGCAGCAACGCCGCGTGAGGGACGACGGCCTTCGGGTTGTAAACCTCTTTTAGTAGGGAAGAAGCGAAAGTGACGGTACCTGCAGAAAAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTATCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAACTGGGGGCTCAACCCCCAGCCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGATTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCAATGGCGAAGGCAGATCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"                    
## [12] "TGGGGAATTTTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGAAGAAGGCCTTCGGGTTGTAAAGCTCTTTTGTCAGGGAAGAAACGGTGTGGGCTAATATCCTGCACTAATGACGGTACCTGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTTGTAAGTCTGTCGTGAAATCCCCGGGCTCAACCTGGGAATTGCGATGGAGACTGCAAGGCTAGAATCTGGCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAACACCGATGGCGAAGGCAGCCCCCTGGGTCAAGATTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [13] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTTAGAGGTGAAAGCCTGGAGCTCAACTCCAGAATTGCCTTTAAGACTGCATCGCTTGAATCCAGGAGAGGTGAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"                         
## [14] "TGGGGAATTTTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGCAGGATGAAGGCCTTCGGGTTGTAAACTGCTTTTGTACGGAACGAAAAAGCTCCTTCTAATACAGGGGGCCCATGACGGTACCGTAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTATGTAAGACAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTTGTGACTGCATGGCTAGAGTACGGTAGAGGGGGATGGAATTCCGCGTGTAGCAGTGAAATGCGTAGATATGCGGAGGAACACCGATGGCGAAGGCAATCCCCTGGACCTGTACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [15] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGTAGGGAACAAGGCCACACGTGGTGTGGTTGAGGGTACTTGCAGAAGAAGCGCCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGCGCAAGCGTTATCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCCGTGAAAGTCCGGGGCTCAACTCCGGATCTGCGGTGGGTACGGGCAGACTAGAGTGATGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCATTAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"         
## [16] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGACGAAGCGAGAGTGACGGTACCTGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTGTCCGGAATTACTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCACGTCGGCTGTGAAAACCCATCGCTCAACGGTGGGCTTGCAGTCGATACGGGCTGACTTGAGTACTGCAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCAGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"                    
## [17] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGGGGGATGAAGGCCTTCGGGTTGTAAACCTCTTTCAGCTCCGACGAATTCAGACGGTAGGAGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCGGCTGTGAAATCCCGTGGCTCAACTGCGGGTCTGCAGTCGATACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"                        
## [18] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGGAGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGAAGAAGTGGTCGGGTTCTCTCGGTCGTGACGGTACCTGCAGAAGAAGCGCCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGCGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAACTCGAGGCTCAACCTCGGGCTTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCCGCAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"      
## [19] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGGGGGATGAAGGCCTTCGGGTTGTAAACCTCTTTCAGCTCCGACGAATTCAGACGGTAGGAGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCGGCTGTGAAATCCCGTGGCTCAACTGCGGGTCTGCAGTCGATACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCAAACA"                        
## [20] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGACGAACACAGACGGTACCTGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTCTGTCGCGTCGGCTGTGAAAACTCGGGGCTCAACCCCGAGCCTGCAGTCGATACGGGCAGACTAGAGTGCTGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCAGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"

5a - Genus-level stacked barchart by DNA_or_cDNA~Fox_Aspect (no filter)

5b - Genus-level stacked barchart by DNA_or_cDNA~Snowpack_Layer (filter)

Oh good. But this need improving. Getting the unclassified component of the families for the top main taxa in it, but without a guarantee for them popping into it. Also, basing it on the >0.1% dataset as per Arwyn’s advice.

ps.final.r.glom.genus.psdf_2 <- data.table(speedyseq::psmelt(ps.final.r.glom.genus))
ps.final.r.glom.genus.psdf_2$Genus <- as.character(ps.final.r.glom.genus.psdf_2$Genus)

ps.final.r.glom.genus.psdf_2[, mean := mean(Abundance, na.rm = TRUE), 
                 by = "Genus"]

ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Sphingomonadaceae"] <- "Unc. Sphingomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Burkholderiaceae"] <- "Unc. Burkholderiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Labraceae"] <- "Unc. Labraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Acetobacteraceae"] <- "Unc. Acetobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Microbacteriaceae"] <- "Unc. Microbacteriaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Deinococcaceae"] <- "Unc. Deinococcaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Gemmatimonaceae"] <- "Unc. Gemmatimonaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Intrasporangiaceae"] <- "Unc. Intrasporangiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Nocardiaceae"] <- "Unc. Nocardiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Ktedonobacteraceae"] <- "Unc. Ktedonobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Cellulomonadaceae"] <- "Unc. Cellulomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          ps.final.r.glom.genus.psdf_2$Family=="Beijerinckiaceae"] <- "Unc. Beijerinckiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &          is.na(ps.final.r.glom.genus.psdf_2$Family)] <- "Unclassified"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) &  is.na(ps.final.r.glom.genus.psdf_2$Family) & ps.final.r.glom.genus.psdf_2$Order=="Frankiales"] <- "Unc. Frankiales"

ps.final.r.glom.genus.psdf_2[(mean <= 0.01), Genus := "Other (<1%)"]

ps.final.r.glom.genus.psdf_2$Genus = factor(ps.final.r.glom.genus.psdf_2$Genus,
                         levels = c("Acidiphilium","Arthrobacter",
                                    "Bacillus",
                                    "Bosea",
                                  "Burkholderia-Caballeronia-Paraburkholderia",
                                    "Cryobacterium","Cupriavidus",
                                    "Delftia","Hymenobacter",
                                    "Labrys","Leptothrix",
                                    "Massilia","Methylobacterium",
                                    "Noviherbaspirillum","Oryzihumus",
                                    "Polymorphobacter",
                                    "Salinibacterium",
                                    "Sphingomonas","Other (<1%)"))
tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477",
               "#4477AA", "#77AADD", "#117777", "#44AAAA", 
               "#77CCCC", "#777711", "#AAAA44", "#DDDD77", 
               "#774411", "#AA7744", "#DDAA77", "#771122", 
               "#AA4455", "#DD7788", "azure3")

ggplot(ps.final.r.glom.genus.psdf_2, 
       aes(x = Month, y = Abundance, fill = Genus)) + 
  facet_grid(DNA_or_cDNA~Snowpack_Layer,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

5c - Genus-level stacked barchart by DNA_or_cDNA~Fox_Aspect (filter)

ggplot(ps.final.r.glom.genus.psdf_2, 
                                   aes(x = Month, y = Abundance, fill = Genus)) + 
  facet_grid(DNA_or_cDNA~Fox_Aspect,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) + 
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

Lest us not forget that size of Other (<1%) corresponds to the different samples included in that category.

I don’t like this, but plotting all samples now, without grouping.

# ps.final.r.glom.genus.psdf_2$Sample_Name <- as.character(ps.final.r.glom.genus.psdf_2$Sample_Name)
# ps.final.r.glom.genus.psdf_2$Sample_Name <- factor(ps.final.r.glom.genus.psdf_2$Sample_Name,
#                                                  levels=unique(ps.final.r.glom.genus.psdf_2$Sample_Name))
ggplot(ps.final.r.glom.genus.psdf_2,
                                   aes(x = Sample_Name, y = Abundance, fill = Genus)) +
  facet_grid(DNA_or_cDNA~Month,
             scales = "free",
             space = "free") +
  geom_bar(stat = "identity",
           position = "fill",
           color="black") +
  scale_fill_manual(values = tol17rainbow,
                    na.value = "gray40") +
  scale_y_continuous(expand = c(0,0),
                     labels = scales::percent,
                     breaks = c(0.25,0.5,0.75,1)) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 35, hjust = 1),
        strip.text.x = element_text(face = "bold",
                                    size = 16),
        strip.background = element_rect(size=1),
        legend.text = element_text(size= 14)) +
  guides(fill = guide_legend(keywidth = 1, keyheight = 1,
                             nrow = )) +
  ylab("Relative Abundance\n")

6 - Unconstrained ordination

nmds_shapes = c(15,16,17,18,6)
nmdscolors_by_month = c("April" = "#A6048B",
                        "June" = "#0704A6",
                        "Early July" = "#727A0B",
                        "Late July" = "#6C0504")
nmdscolors_by_snpl = c("Top" = "#151AC6",
                       "Mid" = "#C6156D",
                       "Basal" = "#890909",
                       "Sup. Ice" = "#091D89",
                       "Glacial Ice" = "#580303")

ps.final.r.dna = subset_samples(ps.final.r, DNA_or_cDNA != "cDNA")
ps.final.r.cdna = subset_samples(ps.final.r, DNA_or_cDNA != "DNA")

ps.final.r.dna.nmds  <- ordinate(ps.final.r.dna,
                         "NMDS",
                         distance="jsd",
                         maxit = 1e3)
## Run 0 stress 0.2107623 
## Run 1 stress 0.2116726 
## Run 2 stress 0.2123577 
## Run 3 stress 0.2133448 
## Run 4 stress 0.2030457 
## ... New best solution
## ... Procrustes: rmse 0.1022482  max resid 0.3790585 
## Run 5 stress 0.2161223 
## Run 6 stress 0.2165051 
## Run 7 stress 0.2139681 
## Run 8 stress 0.2087422 
## Run 9 stress 0.1999096 
## ... New best solution
## ... Procrustes: rmse 0.05625804  max resid 0.2574247 
## Run 10 stress 0.2102544 
## Run 11 stress 0.2150366 
## Run 12 stress 0.2026495 
## Run 13 stress 0.2154225 
## Run 14 stress 0.2213739 
## Run 15 stress 0.2098726 
## Run 16 stress 0.2163894 
## Run 17 stress 0.2031956 
## Run 18 stress 0.2015388 
## Run 19 stress 0.2044157 
## Run 20 stress 0.212047 
## *** No convergence -- monoMDS stopping criteria:
##     20: stress ratio > sratmax
ps.final.r.cdna.nmds  <- ordinate(ps.final.r.cdna,
                         "NMDS",
                         distance="jsd",
                         maxit = 1e3)
## Run 0 stress 0.1654426 
## Run 1 stress 0.1685203 
## Run 2 stress 0.1630985 
## ... New best solution
## ... Procrustes: rmse 0.02948222  max resid 0.1409178 
## Run 3 stress 0.1644532 
## Run 4 stress 0.1734149 
## Run 5 stress 0.1665583 
## Run 6 stress 0.1680772 
## Run 7 stress 0.1626706 
## ... New best solution
## ... Procrustes: rmse 0.120582  max resid 0.4752591 
## Run 8 stress 0.1623674 
## ... New best solution
## ... Procrustes: rmse 0.09400929  max resid 0.3635666 
## Run 9 stress 0.166825 
## Run 10 stress 0.1621902 
## ... New best solution
## ... Procrustes: rmse 0.01076304  max resid 0.0346839 
## Run 11 stress 0.161498 
## ... New best solution
## ... Procrustes: rmse 0.06954823  max resid 0.2343931 
## Run 12 stress 0.1714477 
## Run 13 stress 0.1642582 
## Run 14 stress 0.1663909 
## Run 15 stress 0.1602385 
## ... New best solution
## ... Procrustes: rmse 0.02810959  max resid 0.1054134 
## Run 16 stress 0.1687445 
## Run 17 stress 0.1691279 
## Run 18 stress 0.1670851 
## Run 19 stress 0.1690936 
## Run 20 stress 0.1652288 
## *** No convergence -- monoMDS stopping criteria:
##     17: stress ratio > sratmax
##      3: scale factor of the gradient < sfgrmin
ps.final.r.dna.nmds.plot <- plot_ordination(ps.final.r.dna, ps.final.r.dna.nmds) +
               geom_point(aes(colour=Month, shape=Snowpack_Layer),
                          size=5, alpha = 0.5, stroke=1.5) +
               scale_shape_manual(values=nmds_shapes) +
               scale_colour_manual(values=nmdscolors_by_month) +
               scale_fill_manual(values=nmdscolors_by_month) +
               stat_ellipse(aes(fill=Month),
                          geom = "polygon",
                          type="t",
                          alpha=0.15,
                          linetype = 2,
                          show.legend = FALSE) +
  ggtitle("DNA nMDS") +
  guides(colour = guide_legend(title="Month"))

ps.final.r.cdna.nmds.plot <- plot_ordination(ps.final.r.cdna, ps.final.r.cdna.nmds) +
               geom_point(aes(colour=Month, shape=Snowpack_Layer),
                          size=5, alpha = 0.7, stroke=1.5) +
               scale_shape_manual(values=nmds_shapes) +
               scale_colour_manual(values=nmdscolors_by_month) +
               scale_fill_manual(values=nmdscolors_by_month) +
               stat_ellipse(aes(fill=Month),
                          geom = "polygon",
                          type="t",
                          alpha=0.2,
                          linetype = 2,
                          show.legend = FALSE)+
  ggtitle("cDNA nMDS") +
  guides(colour = guide_legend(title="Month"))
plot_grid(ps.final.r.dna.nmds.plot + theme_bw(), 
          ps.final.r.cdna.nmds.plot + theme_bw())

7 - Alpha-diversity metrics.

month_comparisons = list( c("April", "June"), c("April", "Early July"), 
                          c("April", "Late July"),
                          c("June", "Early July"), c("June", "Late July"),
                          c("Early July", "Late July"))

site_comparisons = list( c("NW_AWS", "N_NE"),
                         c("NW_AWS", "SWS1SE"),
                         c("N_NE", "SWS1SE"))

slayer_comparisons = list( c("TOP", "MID"), c("TOP", "BASAL"), c("TOP", "SUP ICE"), 
                           c("TOP", "GL ICE"),
                           c("MID", "BASAL"), c("MID", "SUP ICE"), 
                           c("MID", "GL ICE"),
                           c("BASAL", "SUP ICE"), c("BASAL", "GL ICE"),
                           c("SUP ICE", "GL ICE"))

All measures for Sample_ID, no stats, facet by DNA/cDNA

plot_richness(ps.final,
              x="Sample_ID", color="Sample_ID",
              measures=c("Observed","Shannon","Simpson")) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Sample_ID") +
  guides(color=guide_legend(ncol=2)) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures by Sample_ID

estimate_richness(ps.final, measures = c("Observed","Shannon","Simpson")) %>% 
  kable(digits = 4) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% row_spec(0,bold=TRUE)
Observed Shannon Simpson
A03.cDNA.1 286 0.7114 0.2103
A03.DNA.1 138 3.2774 0.9407
A04.cDNA.2 163 2.0361 0.7852
A04.DNA.2 75 1.6409 0.6083
A05.cDNA.3 515 2.2067 0.7048
A05.DNA.3 1390 5.1230 0.9813
A06.cDNA.4 125 1.3510 0.6058
A06.DNA.4 210 3.9337 0.9650
A07.cDNA.5 1635 4.6591 0.9488
A07.DNA.5 609 3.3602 0.8726
A08.DNA.6 563 4.0141 0.8786
A09.cDNA.7 1875 5.5902 0.9807
A10.cDNA.8 842 4.7339 0.9520
A10.DNA.8 179 3.3916 0.9493
A11.cDNA.9 3 0.2518 0.1121
A11.DNA.9 1819 5.3121 0.9764
A12.cDNA.10 54 1.4929 0.6133
A12.DNA.10 2719 5.7179 0.9837
B01.cDNA.11 453 4.0612 0.9622
B01.DNA.11 187 3.5277 0.9546
B02.cDNA.12 733 4.5851 0.9683
B02.DNA.12 310 2.1151 0.7154
B03.cDNA.13 38 1.1241 0.4091
B03.DNA.13 457 4.2859 0.9559
B04.cDNA.14 67 1.1867 0.4139
B05.cDNA.15 1307 5.7881 0.9881
B06.cDNA.16 134 1.4982 0.6595
B06.DNA.2 156 0.6586 0.2271
B07.DNA.3 1302 5.3971 0.9910
B08.cDNA.18 21 0.9250 0.5154
B08.DNA.4 986 5.2137 0.9805
B09.DNA.5 105 2.9238 0.9150
B10.cDNA.21 400 4.5700 0.9771
B10.DNA.6 377 4.1547 0.9598
B12.cDNA.23 29 0.9625 0.3511
B12.DNA.8 1998 5.3382 0.9817
C01.DNA.9 165 3.3370 0.9306
C02.cDNA.26 907 4.1753 0.9419
C03.cDNA.27 537 3.3393 0.8840
C03.DNA.B2.1 677 5.1429 0.9887
C04.cDNA.28 307 3.3747 0.8749
C04.DNA.2 715 4.2502 0.9430
C05.DNA.3 1739 5.2347 0.9838
C06.DNA.4 236 3.4402 0.9006
C07.DNA.5 2789 5.9749 0.9894
C08.DNA.6 307 4.5001 0.9837
C09.cDNA.33 36 1.3284 0.5802
C09.DNA.7 148 3.8243 0.9690
C10.DNA.8 383 3.2676 0.8801
C11.cDNA.35 65 1.5747 0.6657
C12.cDNA.36 19 0.4213 0.1462
C12.DNA.10 396 2.2692 0.7492
D01.DNA.B4.1 147 0.6873 0.2422
D02.cDNA.38 128 2.1030 0.7780
D03.cDNA.39 266 3.4007 0.9035
D04.cDNA.40 79 1.5083 0.6114
D05.cDNA.41 154 1.8281 0.7377
D05.DNA.5 90 0.6634 0.2169
D06.cDNA.42 841 3.6042 0.9024
D06.DNA.6 2292 5.5649 0.9810
D07.cDNA.43 1444 4.5026 0.9489
D07.DNA.7 467 1.8746 0.5308
D08.cDNA.44 120 2.6011 0.8428
D08.DNA.8 1045 4.1510 0.9245
D09.DNA.9 2986 5.9104 0.9895

All measures for Month, stats included, facet by DNA/cDNA

plot_richness(ps.final,
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Month") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Month, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed

plot_richness(subset_samples(ps.final, 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)),
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Month") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Month, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)

mini_month_comparisons = list(c("April", "June"), 
                          c("April", "June"))

plot_richness(subset_samples(ps.final, Month != "Early July" & Month != "Late July"),
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Month") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Fox_Aspect, no stats, facet by DNA/cDNA

#now without stats, only shapes
plot_richness(ps.final,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Fox Aspect") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Fox_Aspect, stats included, facet by DNA/cDNA

plot_richness(ps.final,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha=0.5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = site_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed

plot_richness(subset_samples(ps.final, 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha=0.5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = site_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
                             Month != "Early July" & Month != "Late July"),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only

plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July")),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed."

plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July"), 
                             Snowpack_Layer != "GL ICE" & 
                             Sample_Name != "30 July_SWS1SE_Slush_R" &
                             Sample_Name != "30 July_SWS1SE_Slush_D" &
                             !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Early July” only, “GL ICE” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
                               Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Late July” only, “GL ICE”, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & Month == "Late July" &
                               Sample_Name != "30 July_SWS1SE_Slush_R" &
                               Sample_Name != "30 July_SWS1SE_Slush_D" & 
                               !(Sample_Name %in% gl_ice_spls)),
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = mini_month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Fox_Aspect") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

All measures for Snowpack_Layer, no stats, facet by DNA/cDNA

#now without stats, only shapes
plot_richness(ps.final,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Snowpack_Layer, stats included, facet by DNA/cDNA

plot_richness(ps.final,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                     symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA

plot_richness(subset_samples(ps.final, Month = c("Early July", "Late July")),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed

plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July") &
                                       Snowpack_Layer != "GL ICE" & 
                                       !(Sample_Name %in% gl_ice_spls) &
                                       Sample_Name != c("30 July_SWS1SE_Slush_R",
                                                        "30 July_SWS1SE_Slush_D")),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning in `!=.default`(Sample_Name, c("30 July_SWS1SE_Slush_R", "30
## July_SWS1SE_Slush_D")): longer object length is not a multiple of shorter object
## length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

All measures for Snowpack_Layer in “Early July”, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & 
                               Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed

## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

All measures for Snowpack_Layer in Late July, stats included, facet by DNA/cDNA, “GL ICE” removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_R” removed, “Gl ice” samples removed

plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & 
                               Sample_Name != "30 July_SWS1SE_Slush_R" &
                               Sample_Name != "30 July_SWS1SE_Slush_D" &
                               Month == "Late July" & !(Sample_Name %in% gl_ice_spls)),
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin(alpha = .5) +
  geom_jitter(position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  labs(color="Snowpack Layer") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations

# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)

8 - Alpha diversity with rarefaction to the smallest number of reads in any sample

ps.final.rare = rarefy_even_depth(ps.final, 
                                  sample.size = min(sample_sums(ps.final)),
                                  rngseed = TRUE, replace = TRUE, 
                                  trimOTUs = TRUE, verbose = TRUE)
## `set.seed(TRUE)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(TRUE); .Random.seed` for the full vector
## ...
## 25951OTUs were removed because they are no longer 
## present in any sample after random subsampling
## ...
plot_richness(ps.final.rare,
              x="Sample_ID", color="Sample_ID",
              measures=c("Observed","Shannon","Simpson")) +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Sample_ID") +
  guides(color=guide_legend(ncol=2)) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Month", color="Month",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = month_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Month") +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(shape=16, position=position_jitter(0.2)) + 
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = site_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Fox Aspect") +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

#now without stats, only shapes
plot_richness(ps.final.rare,
              x="Fox_Aspect", color="Fox_Aspect",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Fox Aspect") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(position=position_jitter(0.2)) + 
  # stat_compare_means(comparisons = my_comparisons)+
  stat_compare_means(comparisons = slayer_comparisons,
                     method = "t.test",
                     var.equal="welch",
                     label.x.npc = "right", 
                     label.y.npc = "top",
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1), 
                       symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

#now without stats, only shapes
plot_richness(ps.final.rare,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Month), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))

plot_richness(ps.final.rare,
              x="Snowpack_Layer", color="Snowpack_Layer",
              measures=c("Observed","Shannon","Simpson")) +
  geom_violin() +
  geom_jitter(aes(shape=Fox_Aspect), size=5, alpha=0.7,
              position=position_jitter(0.2)) + 
  facet_wrap(DNA_or_cDNA~variable, scales = "free") +
  labs(color="Snowpack Layer") +
  scale_shape_manual(values=nmds_shapes) +
  theme(axis.text.x = element_text(angle = 20,
                                   vjust = 1,
                                   hjust = 1,
                                   size = 12),
        axis.title = element_blank(),
        strip.text = element_text(size = 12, margin = margin()))